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We propose a new Monte Carlo scheme to study the late- 
time dynamics of a 2-dim hard sphere fluid, modeled by a 
tethered network of hard spheres. Fluidity is simulated by 
breaking and reattaching the flexible tethers. We study the 
diffusion of a tagged particle, and show that the velocity auto- 
correlation function has a long-time t~ x tail. We investigate 
the dynamics of phase separation of a binary fluid at late 
times, and show that the domain size R(t) grows as t 1 ^ 2 for 
high viscosity fluids with a crossover to t 2 ^ 3 for low viscos- 
ity fluids. Our scheme can accomodate particles interacting 
with a pair potential V(r), and modified to study dynamics 
of fluids in three dimensions. 

PACS: 64.60.Cn, 64.60.My, 64.70.-p, 61.20.Ja, 51.10.+y 

With the advent of modern computers, there have been 
several large scale molecular dynamics (MD) |IJ , lattice- 
gas automata (LG) and Langevin simulations (LS) || 
of the late-time dynamics of fluids. These simulations 
have been used to compute dynamical correlation func- 
tions and to study fluid flow in different geometries. MD 
and LS have also been used to study the dynamics of 
phase ordering of binary fluids, where the relative con- 
centration of the fluids is coupled to the fluid velocity. In 
many such applications, one is interested in late-time hy- 
drodynamical behaviour. Since MD provides an accurate 
description of the microscopic physics, it may not probe 
dynamics at such late time, unless one makes a sizeable 
computational investment. LS are a more coarse-grained 
description ; however the dynamical equations are too 
complicated to solve, and it is not clear that the vari- 
ous simplifying approximations made, do not miss out 
features which might affect the physics at large length 
and time scales. 

There is a considerable advantage if Monte Carlo (MC) 
simulations could be used instead. MC simulations Q 
have been remarkably successful in the study of the dy- 
namics of alloys, magnets and so on, but have never been 
used in dynamical studies involving fluids. This is be- 
cause there is no natural way to incorporate the momen- 
tum density in a MC simulation. It would be desirable to 
devise an MC scheme to handle fluidity, since it is more 
coarse-grained than MD and more microscopic than LS. 
More recently Lattice Boltzmann (LB) || simulations 
have been employed with this specific goal in mind, but 
since such simulations do not depend on the specific form 
of the hamiltonian, one has to choose the equilibrium dis- 
tribution functions consistent with interfacial profiles and 



equilibrium densities. 

For the stated reasons and more, we present in this 
Letter, a novel MC technique which successfully simu- 
lates the fluidity of a fluid. Though the algorithm can 
be extended to dimensions d > 2, and to any pair poten- 
tial, we demonstrate its efficacy for a hard sphere fluid 
in d = 2. We determine the single-particle diffusion coef- 
ficient and the long-time tail of the velocity autocorrela- 
tion function for dense fluids. We also study the kinetics 
of phase ordering of a binary fluid in 2-d and compute the 
concentration correlation function to extract the growth 
laws described below. 

Our model 2-d fluid consists of hard spherical beads 
(vertices) (the total number of beads AT, is fixed), re- 
stricted onto a finite region of area A oilZ 2 . All the beads 
are linked together by straight flexible tethers (bonds), 
in such a way as to triangulate this region. The teth- 
ers do not intersect each other — each configuration is 
thus a connected planar graph. We ensure that the lo- 
cal coordination number of every bead lies between 3 
and 9. Since the particles are restricted to a 2-d plane 
with zero curvature, the distribution of local coordination 
numbers should be symmetric about 6. The beads are in- 
finitely repulsive at distances less than a bead diameter 
Imin = a (for all vertex pairs) or greater than l max (for 
vertices connected by tethers), and so the local tether 
length can vary between these two limits. Our choice of 
lmax = 5v3a, guarantees that half of the attempted MC 
moves are accepted. We distinguish between external 
(edge) vertices (V E ) and internal (bulk) vertices (V 1 ). 
An external (edge) tether (T E ) connects two external 
vertices, while an internal (bulk) tether (T 1 ) connects at 
least one internal vertex. 

Our Monte Carlo simulation should allow for config- 
uration changes which sample all permissable regions in 
phase space. The configuration changes of our 2-d fluid 
consist of the movement of beads with fixed connectiv- 
ity (bead moves) and the movement of tethers with fixed 
vertex positions (flip moves). As we shall see below, the 
flip moves are crucial in maintaining fluidity. We shall 
use rigid boundary conditions (on a hexagonal frame), 
and so neither V E nor T E are moved. 

The bead moves are effected by randomly choosing a 
bead i and then translating i to a random point (with a 
uniform distribution) within a square (of size I = 2a) cen- 
tered on the old position of i. The movement is accepted 
if the new bond lengths lie between l m i n and l ma x and 
the graph remains planar (i.e., no bond intersections). 
With each particle i, we can identify a unique n-gon 
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(concave or convex) whose vertices {v±, V2, ■ ■ ■ , v n }, ar- 
ranged in cyclic order, are connected to i. Let us de- 
note the areas of the triangles {iv%V2, ■ ■ ■ , iv n -xv n } as 
{A{iv\V2), ■ ■ ■ , A(iv n —iv n )}. Planarity is maintained if 
A(iv\V2) + . . . + A(iv n —iv n ) remains unchanged after the 
move on i. Our unit of time is set by one Monte Carlo 
sweep (MCS), defined as N attempted bead moves. 

Fluidity can be generated by a combination of bead 
and flip moves. The flip moves are a modification of 
the bond reconnection algorithm, developed to study the 
equilibrium behaviour of fluid membranes || . A bond t^ 
(connecting vertices i and j) is picked at random. With 
every internal bond ty , we can identify two triangles ijv\ 
and ijv2 on either side of tij which have the smallest 
area. This defines a quadrilateral, with i and j being a 
pair of opposite vertices. The bond tij is now flipped, 
only if v\ and V2 are not already connected by a bond, 
so that it now connects v\ and V2 (leaving the vertices 
i and j unconnected by a bond). This flip is accepted 
provided the length of t VlV2 is less than the maximum 
allowed length l max , and if t vlV2 does not intersect any 
other tether (planarity). Note that the total number of 
bonds is conserved during this operation. During one 
MCS, we make Nfu p attempts at flipping bonds. 

How well does this algorithm mimic the statics and 
dynamics of a simple fluid ? The equilibrium proper- 
ties of a simple fluid can be derived from the parti- 
tion function Z = Tr s Tr p exp(—f3F[{g},{p}]) (where 
F = J [g 2 /2po + V({p})} is the free-energy functional 
of the local momentum density g(r, t) and the mass den- 
sity p(r, t) of the fluid). In our simulation, we define 
the local velocity for the i-th particle, = gi/m, as 
its vector displacement in a single MCS. Even though 
particle displacements are picked out with a uniform dis- 
tribution from a square of side I centred at the original 
position of the particle, large displacements, which have 
a high chance of violating the tethering and planarity 
constraints, are suppressed. Clearly, the larger the co- 
ordination number of the i-th particle, the higher is the 
chance that large displacements are suppressed and so 
the smaller is The particle displacements across the 
whole sample are independent of each other and so by the 
central limit theorem, the displacements (and hence the 
velocities) of particles tends to a normal (Maxwellian) 
distribution over several MC times. The time scale over 
which this happens is the velocity equilibration time r u . 
From equipartition, a computation of {ui 2 ) gives the ki- 
netic temperature 2fceT/m where m is the mass of the 
particle. On the other hand, the trace over p is identical 
to a sum over configurations with (dynamically) varying 
local coordination number. Since the tethering is dynam- 
ical, there is no restriction on the allowable configurations 
sampled by our fluid. 

Equilibrium is achieved through collisions between par- 
ticles which result in velocity exchanges. In our simula- 
tions, "collisions" occur either through the rejection of 
certain MC moves due to the hard sphere constraint, or 
through flip moves. Collisions resulting from flip moves, 



transform particles with low coordination number to high 
and vice versa. From our discussion above, this results 
in an exchange of velocities. It is clear that particle mo- 
menta (and angular momenta) are not conserved during 
an elementary "collision" , and so our algorithm will not 
mimic fluid dynamics at short time scales. However, av- 
eraged over several collisions, both momenta and angu- 
lar momenta can be seen to be conserved. Just as in 
a Langevin simulation, momentum is conserved only in 
a statistical sense. Our MC dynamics thus mimics the 
dynamics of a fluid over several "collision" time scales. 
Starting with a random triangulation (which we gener- 
ate by repeated bead and flip moves on an initial trian- 
gular lattice), we allow the fluid to evolve via our MC 
algorithm. We measure various dynamical transport co- 
efficients after discarding all configurations arising from 
typically the first 10 Monte Carlo sweeps. The single- 
particle diffusion coefficient D s can be extracted by mea- 
suring the mean-square displacement of a tagged particle 
(averaged over several initial conditions and particles) for 
different values of Nju p and density (Fig. 1). Not sur- 
prisingly, D s (and hence the inverse viscosity t]^ 1 ) is an 
increasing function of Nfu p (inset Fig. 1) at constant 
density, before it saturates. 
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FIG. 1. The single particle diffusion coefficient D s at 
different Nfu p , for Nfu p = TV, 2N, 5N and 97V going from 
bottom to top, is given by the value of (r 2 )/4i at the plateau. 
The inset shows D s as a function of Nfu p (in units of TV). 
The area fraction Niva 2 /A = 0.14, is kept fixed. 

A crucial test of whether our MC algorithm correctly 
describes the hydrodynamic behaviour of fluids, is the 
occurence of "long-time tails" in the measured velocity 
autocorrelation function Z(t) =< Uj(0) • Uj(£) > of dense 
(or highly viscous) fluids 0. We compute Z{t) (averaged 
over several initial conditions and particles) for different 
values of Nfu p and density. We find a convincing i _1 
tail (Fig. 2) in the Z{t) of high density (large N) and 
high viscosity (low N*n p = 0.25./V) fluids. The inset of 
Fig. 2 shows the exponential decay (Enskog result) of 
Z(t) for lower viscosity (Nm p = 10N). Note that the 
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early time fall from Z(Q) — 1, does not show up, since 
for the densities under consideration, several "collisions" 
have taken place within one MC time step. The inset 
also displays the t~ x tail for comparison. 
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FIG. 2. The velocity autocorrelation function Z(t) for 
high density, high viscosity fluids exhibits a clear t^ 1 long- 
time tail. Inset : Z(t) decays to zero for low viscosity fluids 
(the i _1 tail is displayed for comparison). 

Encouraged by this success, we study the dynamics 
of phase separation of a binary fluid using our MC al- 
gorithm. Let us recount that when two immiscible flu- 
ids, like water+toluene, are cooled below their coexis- 
tence curve, they phase segregate into water-rich and 
toluene-rich regions, separated by sharp interfaces. At 
late times, the system enters a dynamical scaling regime 
H , where the equal-time concentration correlation func- 
tion behaves as g(r, t) = g(r/t z ). The growth exponent z 
is independent of microscopic details and depends on the 
existence of conservation laws. The scaling form dchncs 
a characteristic length scale R(t) ~ t z , which measures 
the typical distance between interfaces. 

The phase separation dynamics of a 2-d binary fluid is 
described by a (conserved) concentration density <p(r, t) 
coupled to a (conserved) momentum density w(r,t). A 
variety of theoretical (based on dimensional analysis 
g,|) and numerical (LS Q and MD @) techniques have 
been used to understand the late stage dynamics of a 2-d 
binary fluid, but have unfortunately given rise to con- 
flicting results. The theoretical analysis and LS, contend 
that, as in 3-d, the growth crosses over from a viscosity- 
dominated, R ~ t, to an inertia-dominated, R ~ t 2 / 3, . 
On the other hand, extensive MD simulations || report 
a late time R ~ t x l 2 growth. 

Our model 2-d binary fluid consists of two types of 
hard spheres, A and B (the total number of A (Na) and 
B (N B ) beads are fixed, N = N A + N B ). The local con- 
centration (j>i , defined as the difference in the densities of 
A and B, takes values +1 (for A) and —1 (for B). The 
dynamics of <f> is described by the usual Kawasaki ex- 
change, which conserves Yli&i- Thus locally <j>i evolves 



by exchanging particles at vertices i and j, where j is 
connected to j by a tether, with a transition probabil- 
ity W(i <-> j) = [1 - taah(AE /2k B T)]/2, where AE is 
the energy difference between the final and the initial 
configuration. The energy is calculated using the Ising 
hamiltonian 



<ij> 



(1) 



where the sum over < ij > is over vertex pairs connected 
by a tether. 

In addition to the two MC moves described above, each 
Monte Carlo sweep now includes N ex attempts at per- 
forming Kawasaki exchanges. Note that the Boltzmann 
weights associated with particle movements and the bond 
flips do not involve H ex . It is clear that fluctuations in 
the local exchange- energy are related to fluctuations in 
the local coordination number, which in turn are related 
to fluctuations in the local velocity. In this way, the con- 
centration variable 4> gets coupled to the velocity u. 
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FIG. 3. Scaling functions for the correlation function 
g(r,t) at (a) high viscosity (Nfu p — 30N) and (b) low viscos- 
ity (Nfup = lOOiV). Here x = r/R(t) is the scaling variable. 

At time t = 0, our 2-d binary fluid is a homogeneous 
(50-50) mixture of A and B in equilibrium at a high tem- 
perature T. We quench to below the critical temperature, 
T < T c w 2.25 J (for a triangular lattice). The homoge- 
neous state is unstable at this temperature, and evolves 
into a final phase separated state by the slow annealing of 
interfaces, conserving the order parameter <f> during the 
process. A time sequence of typical configurations of <j> 
(resembling those in Ref. ||), clearly exhibit statistical 
self-similarity at late times. As a quantitative measure, 
we compute the circularly averaged pair correlation func- 
tion g(r, t) = N~ x X) x (^( x + r ' O0( x > £))) averaged over 
several realisations of the initial configurations of <f> (it 
was sufficient to average over 10 realisations for good 
statistics). The domain size R(t) is extracted from the 
first zero of the correlation function, g(R(t), t) — 0. The 
correlation function, g(r,t) exhibits dynamical scaling, 
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g(r,t) — g(r/R(t)), at late times (Fig. 3), both for high 
viscosity Nm p = SON (Fig. 3a) and for low viscosity 
NfUp = 1007V (Fig. 3b) fluids. Our study seems to indi- 
cate that the scaling function depends on the viscosity of 
the fluid. 




FIG. 4. Domain size (log-log plot) R(t) scales as (o) t 1 ^ 2 
for high viscosity (Nfu p = 307V) and (+) t 2 ^ 3 for low viscosity 
(TV/Hp = 100JV). The lines with slopes 1/2 and 2/3 are aids 
to the eye. 

The domain size R(t), determined both from the corre- 
lation function and the interfacial energy density, is found 
to scale as R(t) ~ t z . As a preliminary check, we perform 
the usual Kawasaki dynamics in the absence of any bond 
flip and obtain a clean t 1 ^ 3 Lifshitz-Slyozov growth. On 
introducing the flip moves that mimic fluidity, we see a 
dramatic crossover to a growth influenced by hydrody- 
namics. For high viscosity fluids, Nfu p — 30N, Fig. 4 
shows a growth consistent with R ~ t 1 / 2 at late times, 
before finite size effects become apparent. A log- log plot 
for N = 7500, shows a z = 0.48 ± 0.03. Simulations for 
larger system sizes and longer times, give the same value 
of z, indicating that this value of the exponent is robust 
upto the late times we have investigated. This growth 
law is in agreement with MD simulations II. On the 
other hand, for low viscosity fluids Nfu p = 100 AT, Fig. 4 
shows a late time exponent z = 0.6 ± 0.05, indicating a 
crossover from the viscosity dominated t 1 / 2 growth to the 
inertia dominated t 2 ! 3 growth. Thus we expect R ~ i 1 / 2 
when f « f x , and R ~ t 2 ! 3 when f > t X) where the 
crossover time t y ~ rj a . Indeed, recent LB simulations 
JlH , have seen the t 2 ! 3 growth for low viscosity fluids. 
At higher viscosities they see an early time i 1//3 growth, 
and claim a crossover to a later time t 2 ! 3 . However their 
published data fL2f clearly indicates a crossover to t 1 ! 2 . 
We claim that the LB results are consistent with our MC 
simulations. As is clear from the above discussion, our 
computational costs are low. This crossover from z = 1/2 
for high viscosities, to z — 2/3 for low viscosities, can be 



understood theoretically, and we defer a discussion to 
another paper |ll[] . 

In this Letter, we have succesfully implemented a 
Monte Carlo algorithm, which describes the late-time 
dynamics of fluids in two dimensions. This algorithm 
incorporates the momentum density of the fluids. We 
have tested this code on a hard sphere fluid and have 
computed the long-time tail of the velocity autocorre- 
lation function. We have studied the phase separation 
dynamics of a binary fluid, and find that the domain 
size R(t) ~ t 1 ! 2 for high viscosity, and crosses over to 
R(t) ~ t 2 ! 3 for low viscosity fluids, with a crossover time 
depending on the viscosity. This algorithm can be easily 
extended to allow for an arbitrary pair potential V{r). In 
a subsequent publication we shall develop an analogous 
MC algorithm to study the late-time hydrodynamics of 
fluids in three dimensions. 

We have enjoyed discussions with Surajit Sengupta, 
Sriram Ramaswamy and Rahul Pandit. 
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